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ABSTRACT 


Prediction  of  demand  and  setting  reorder  levels  by  ex¬ 
ponential  smoothing  and  mean  absolute  deviation  is  studied 
for  an  inventory  system.  The  relation  between  the  stand¬ 
ard  deviation  and  mean  absolute  deviation  of  demand  dis¬ 
tribution  is  discussed.  Asymptotic  first  and  second 
moments  of  the  exponentially-smoothed  forecast  errors  (MAD) 
are  derived.  Results  obtained  from  simulation  of  several 
normal  systems  indicated  that  the  smoothing  technique  is 
inferior  to  the  classical  maximum  likelihood  estimation 
me thod. 
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I .  INTRODUCTION 


In  recent  years  the  technique  of  forecasting  time 
series  by  exponential  smoothing  has  received  a  considerable 
amount  of  attention  in  the  literature.  As  a  forecasting 
technique  in  inventory  theory,  it  is  particularly  suited  to 
the  situation  wliere  demands  are  time  dependent  and  corre¬ 
lated  for  much  of  its  intuitive  appeal  rests  on  weighting 
the  past  in  such  a  way  that  the  most  distant  past  receives 
less  and  less  weight.  VHien  demands  are  not  so  related, 
smoothing  may  or  may  not  be  a  suitable  forecasting  tech¬ 
nique.  In  particular,  when  demands  are  independent  and 
from  a  common  distribution  such  as  for  the  so-called  fast 
moving  items,  more  classical  methods  like  maximum  likeli¬ 
hood  may  be  preferable  because  of  the  optimal  properties 
that  are  often  demonstrable.  Even  in  such  circumstances, 
however,  exponential  smoothing  is  still  a  candidate  for  a 
forecasting  scheme  and  it  may  very  well  be  tliat,  while  not 
optimal  ,  the  difference  in  such  criteria  as  stockout  risk 
for  example,  may  be  negligible  compared  to  optimal  proce¬ 
dures  from  a  practical  point  of  view.  In  that  case,  fore¬ 
casting  rules  need  not  be  changed  for  items  differing  in 
demand  pattern  and  a  certain  uniformity  may  thus  be 
achieved. 

The  main  purpose  of  this  thesis  is  to  investigate  a 
model  of  an  inventory  system  used  by  the  Naval  Supply 
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System  Command  (NAV  SUP)  following  Brown  [Ref.  l].  The 
model,  as  described  by  Ref.  1  and  Ref.  5,  is  briefly  as 
follows : 

Single  exponential  smoothing  is  an  operator  given  in 
Ref.  1  as 

t-1 

s,(x)  =.  a  I  U’'  (1-1) 

k=0 

where  is  the  total  demand  of  the  i^^  period 

0(  is  the  smoothing  constant  with  0  <  a  <  1 

=  1  -  a 

Intuitively  speaking  the  smoothing  operator  is  a 
method  of  forecasting  by  weighting  past  demands  in  such  a 
way  that  the  most  recent  demand  is  weighted  by  a  and  the 
remaining  demands  are  given  less  and  less  weight  as  time 
progresses. 

The  standard  deviation,  (T,  of  the  demand  distribution, 
which  is  needed  to  set  the  reorder  levels,  is  estimated  by 
the  following  three  statistics; 

e^  =  -  Sj._^(x),  t*"^  period  forecast  error, 

t^l 

A=  a  \  !  ej  fl-2) 

k=0 
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estimate  of  the  mean  absolute  deviation  of  forecast  errors, 
and 


estimate  of  the  standard  deviation  of  the  demand  distribu¬ 
tion. 

This  estimate  of  a,  as  explained  in  Ref.  1,  is  based 
on  the  mean  absolute  deviation  (MAD) ,  A,  of  the  forecast 
errors.  Reorder  levels,  based  on  the  predetermined  risk 
level  j|)  (defined  to  be  the  probability  of  running  out  of 
stock),  are  found  as  follows: 

ji)  =  P(X  >  4  +  ka) 

Here,  4  is  the  mean  of  X,  and  k  is  a  constant  which 
can  be  determined  from  the  probability  tables  for  a  given 
risk  level,  ij).  The  reorder  level  is  then  given  by 

Reorder  Level  =  4  +  kCT  (1-3) 

where  k  is  the  given  above. 

In  the  analysis  this  model  was  slightly  modified  be¬ 
cause  of  the  following  considerations: 

1.  X^,  the  demand  asked  from  the  inventory  system 

before  it  is  placed  in  operation,  is  zero.  Although  one 
may  predict  a  number  for  a  fictitious  inventory,  there  will 
be  no  observation  for  X  .  Hence  e  is  zero. 
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2.  0^  and  can  be  assvuned  to  be  approximately  equal 

to  zero  for  values  of  N  >  30,  if  0.01  s  a  s  0.99.  In  view 
of  this,  the  modified  form  of  (1-1)  and  (1-2)  used  in  the 
analysis  was 

t-1 

S^(x)  =  a  ^t-k 

k=0 

t-1 

i  =  a  ^  ^  1  ej._^|  (1-5) 

k=0 

In  this  model,  not  only  is  demand  forecasted  by  expo¬ 
nential  smoothing  but  demand  variability  is  estimated  by 
means  of  exponentially  smoothed  absolute  deviations.  Such 
estimates  are  in  turn  used  to  set  safety  levels  for  inven¬ 
tory  order  rules. 

The  model  described  above  has  been  subjected  to  many 
criticisms.  Zehna  [Ref.  7 3  has  indicated  that  the  estimate 
of  standard  deviation  might  be  very  poor  depending  upon  the 
underlying  demand  distribution.  Based  on  the  results  ob¬ 
tained  by  Zehna  [Ref.  Sj,  the  analysis  in  this  paper  takes 
on  its  starting  point  on  investigation  of  the  ratio  A/a, 
a  key  element  in  the  estimation  of  ff.  It  is  quite  evident 
that  any  weakness  in  the  estimation  procedure  for  a  affects 
the  reorder  levels.  An  underestimate  imposes  higher  risk 
and  an  overestimate  causes  unnecessary  expenditures.  In 
particular,  the  resulting  consequences  could  be  catastroph¬ 
ic  in  a  military  inventory  system. 
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It  is  quite  evident  that  the  probability  distributions 
of  certain  key  statistics  would  be  desirable  in  any  exten- 
tion  of  the  existing  analysis.  Although  some  of  the  dis¬ 
tributions  required  have  been  derived,  the  distribution  of 
t 

A  ,  a  most  important  statistic,  could  not  be  derived  ex¬ 
plicitly.  However,  the  first  and  second  asymptotic  moments 
of  A  were  obtained,  and  it  is  shown  that  the  asymptotic 
variance  of  A  is  finite  for  the  case  where  the  demand  has 
a  normal  distribution.  Therefore,  the  next  step  of  the 
analysis  used  simulation  techniques  to  gain  an  insight  into 
the  nehavior  of  the  distribution  of  A  ,  and  to  observe  the 
effects  of  a  finite  variance.  Reorder  levels  and  their 
dispersions  along  with  the  risk  levels  attained  are  analyz¬ 
ed,  The  resulting  errors  due  to  the  variability  of  safety 
levels  is  examined  in  detail  and  the  consequences  are 
evaluated  from  several  points  of  view  when  demands  are 
normal.  Moreover  the  results  are  compared  to  maximum  like¬ 
lihood  methods  establishing  that  a  considerable  difference 
exists  in  the  two  methods.  For  items  satisfying  the  demand 
assumptions  considered  here,  then,  it  is  established  quan¬ 
titatively  that  exponential  smoothing  is  definitely  inferior 
to  more  classical  methods. 
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II.  RATIO  OF  A/g 


The  mean  absolute  deviation  of  a  probability  distri¬ 
bution  is  given  by 

A  =  E(|  X  -  mI  )  . 

For  a  normal  distribution  the  ratio  A/cr  is  the  con¬ 
stant  =0.8.  Brown  [Ref.  l]  in  his  work  on  smoothing 

and  prediction  states, 

"Hence,  the  mean  absolute  deviation  is  proportional  to 
the  standard  deviation,  and  the  ratio  may  depend  on  the 
form  of  the  distribution  p(t),  but  dependence  is  slight." 

Now,  as  indicated  by  Zehna  in  [Ref.  7],  this  is  not 
always  the  case.  For  example,  consider  a  Poisson  distri¬ 
bution  with  parameter  X.  A/a  is  then  a  function  of  X  and 
it  may  be  much  smaller  than  0.8  for  small  values  of  X 
[Ref. 8].  However,  there  are  some  distributions  where  the 
approximation  0.8  for  the  ratio  A/a  is  a  good  one  under 
certain  conditions.  The  quote  above  is  a  rather  general 
statement  which  may  be  justified  as  follows; 

Let  X  be  a  continuous  random  variable  with  probability 
density  function  (p.d.f.)  f(x)  and  let  u  and  _1  be  the  upper 
and  lower  limits,  respectively,  of  the  random  variable  X  on 
the  real  line.  Then 

u 

A  =  E{!x-)i|)  =  flx-ul  f(x)dx 

^  I 
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T 


7 


7 


or 


M 


dx 


Hence, 


A  =  2 


F(x)dx 

t 


Here,  F(x)  is  the  cumulative  distribution  of 

Then 


A 

cr 


Let  X  =  ay  +  n.  Then, 


A 

a 


u 

2  I  F(oy  +u)dy 


Now  suppose  that  F(x)  is  any  distribution  function  such 
that 


F(x) 


00 


where  $(z)  is  the  cumulative  standard  normal  distribution 
function.  Then  for  large  values  of  n/a 


F(ay  +  fi)  =  <I>(y) 
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III.  PROBABILITY  DISTRIBUTION  OF  A 


AND  OTHER  ASSOCIATED  STATISTICS 

Of  central  interest  in  the  use  of  MAD  is  the  ability 

to  estimate  the  demand  variability.  The  variance  of  the 

forecast  error  is  always  a  function  of  the  variance  of  de- 
2 

mand,  a  .  Therefore  the  MAD  of  the  forecast  errors  is  also 
a  function  of  a.  Since  MAD  is  a  linear  function  of  the 
standard  deviation  [Ref.  l]  and  if  MAD  of  the  forecast 
errors  is  estimated,  then  a  can  be  estimated  as  described 
in  [Ref.  l].  But  in  this  estimation  there  are  two  points 
to  be  considered. 

1.  Invariance  Principle; 

a  and  A  are  related  to  each  other  by  a  constant, 

X  •  6  •  f 

a  =  cA  . 

But  this  relation  is  in  the  population,  and  it  does  not 
mean  that  the  estimator  of  a  and  A  will  enjoy  the  same 
relationship  necessarily.  If  these  estimates  mentioned 
above  were  MLE,  then  the  invariance  principle  does  hold. 
However,  the  definition  of  A  indicates  that  this  is  not 
the  case. 

2 .  The  Constant  of  Proportionality  of  A/a: 

If  the  underlying  demand  distribution  is  not 
normal,  then  the  assumption  that  the  proportionality  con¬ 
stant  is  approximately  0.8  needs  to  be  validated.  For  this 
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case  the  current  model  can  be  applied  by  a  slight  modifica¬ 
tion  on  a,  that  is,  replacing  by  the  reciprocal  of 

the  appropriate  proportionality  constant. 

2 

In  this  section  we  assume  that  demand  is  Normal  (a, a  ), 

where  a  is  a  constant.  Derivation  of  the  probability  dis- 
■i 

tribution  of  A  requires  probability  distribution  for  the 
following  statistics: 

A.  EXPONENTIAL  SMOOTHING  OF  PAST  OBSERVATIONS 

(Estimate  of  the  Expected  Value  of  Demand) 

t-1 

s^(x)  =  a  ^ 

k=0 


Since  all  X's  are  mutually  independent,  the 

distribution  of  S^(x)  for  large  values  of  t  is  known  to  be 

2  OL 

a  normal  distribution  with  mean  a  and  variance  a  (1  + 


B. 


FORECAST  ERRORS 


t-i-2 

>  .  =  x^  .  -  a  )  x^  .  ,  , 

t-1  t-1  L,  t-i-k-1 

k=0 


i  <  t 


e^_j^  has  a  normal  distribution  with  the  follow¬ 
ing  mean  and  variance ; 

4=0  for  large  values  of  t; 


V(et_i)  = 


(1  + 


l+)9 


for  large  values  of  t. 
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C.  JOINT  DISTRIBUTION  OF  FORECAST  ERRORS 

Let  U  and  V  be  t-dimensional  vectors  defined 


as 


—  — 

—  — 

e. 

u. 

t 

t 

®t-l 

• 

^t-1 

• 

• 

• 

• 

• 

• 

U  = 

• 

• 

• 

• 

• 

®2 

•  •  <N 

"l  1 

where  ...»  are  constants.  Let  Y  be  defined  as 

T 

Y  =  U  V 


or 


t 


i=l 


Substituting  e^  in  terms  of  the  x's, 
t  i=2 

i=l  k=0 


=  I  “A  -  “  I  i  ='i-k-i 

i=l  i=l  k=0 
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Hence, 


t-1  i-1 

=  Z  C“t-i  -  “  Z  “t-k ) 

i=0  k=-0 

Since  the  coefficients  of  the  X‘s  are  constants  and  the  X's 
are  mutually  independent  normal  variables,  it  follows  that 
Y  has  a  normal  distribution.  Hence  by  the  definition  of  a 
multivariate  normal  distribution  [Ref.  2],  the  forecast 
errors  have  a  multivariate  normal  distribution.  The  corre¬ 
lation  coefficient  between  e.  .  and  e^  .,  where  i  <  j,  can 

t-i  t-j 

be  computed  as  follows. 

Let 

Y  =  ®t-j  where  i  <  j  <  t. 


V(Y)  =  V(ej._^)  +  (V(ej._^)  +  2Cov(e^._^,  e|._ .)  , 


or 


Cov(ej._.  ,e^._  .)  =  I  lv(Y)  -  V(e^_  . ) -V  (e^_^ )  J  (3-1) 


Rewriting  in  terms  of  the  X's, 


Y  =  X.  .  -  ax  -  a/9X  ,  - 

t-i  t-i-1  t-i-2 


t-j 


^x 


t-j-1 


+  X  .  -  aX  .  _  aj9X 

t-D  t-D-1  t-j-2 


20 


Rearranging  terras  on  the  right-hand  side 


If  =  ='t-i  -  “  I  ^  - 

k=0 

t-^j-1 

k=0 


Since  the  terms  in  the  above  expression  are  rautually  inde¬ 
pendent,  the  variance  of  Y  for  large  values  of  t  can  be 
found  as  follows: 


V(Y)  =  0“^  (  1  +  a 


0 


1  fl2(j-i-l)  . 

^  "v  +  (l-a/SP  ^  ^ 

i-r 


2 

)  + 


Substituting  V(Y),  into  (3-1), 


Cove(ej._^,e^_j)  =xC^‘‘'l+|  d-l*^  ^  )  + 


(l-a/9 


j-i-1, 


l+/i  y 


Hence,  the  asymptotic  correlation  coefficient  between  e^  ^ 

and  e,  .  becomes 
t-D 


(l-a?3  1 


(3-2) 
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Now  consider  the  following  function. 


Here  y  is  a  continuous  variable.  Then 
do*  a  f  2B^ 

dy  “  4  V"  ^2  ^  +  2^  log  p  +  2$^^  log  $  J  ^ 

l±Af'_  2a^y  ^ 

4  V  j9  log  ^  +  - 1 —  log  Pj 


Rearranging  the  terms. 


a£  .  Jl+^)a  ,  a2(i+s),y  s 
2  2?  +  2,2  ■  ; 


f '  =  (a  -  -ii^)  =  -  loc.  H 

Now  log  />  is  negative:  therefore,  the  first  derivative  of 
P’,  with  respect  to  y  ie  positive.  Hence  p',  is  a  monotoni 
increasing  function  of  y.  Then  positive  integral  values  of 
y  also  lie  on  this  same  continuous  curve.  Actually,  p  co¬ 
incides  with  p',  whenever  y  is  a  positive  integer.  Hence  p 
is  a  monotonic  increasing  function  of  j-i.  The  supremum  of 
P  is  found  by  letting  j-i  in  (3-2) ; 

sup  p  =  (_1  +  ^  (111)  iij.  1 
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The  minimum  value  of  p  is  found  by  letting  j-i  =  1: 
*>min  =  ^  Cl*  i?|  (d-l)  *  (1*^)^)  *  (l-“)^)  -  1 

=  Cl  +  (1+^)®  +  0^^  -  1 

=  i^- 1 


The  minimum  value  of  $  may  be  zero;  hence,  p  .  =  -0.5  . 

min 

D.  DISTRIBUTION  OF  1  .| 

t-i 

Rewriting  |  e^_^|  as 


The  term 


t-i 


-''<1*  if?> 


has,  asymptotically,  the  chi-square  distribution  with  one 
degree  of  freedom  [Ref.  3]. 


1 


The  term 


has,  asymptotically,  the  chi  distribution.  Hence  the  mean 
value  and  variance  of  is  given  by 


'(  I  e^._^  I  )  ^/  n 


(3-3) 


V(le,.,l  )  (1.^, 


(3-4) 


An  attempt  to  derive  the  probability  distribution  of  was 
made  utilizing  the  change-of-variable  technique,  and  by 
defining  the  following  transformation: 

t-1 


A  =  a  I  I 


k=0 


^1  ^  ®t-l 


^t-r  «i 


The  inverse  transformation  becomes 


e  =  ±  Y 
®1  ^t-1 


®2  ^t-2 


t-1 


k=l 

This  transformation  is  obviously  not  one-to-one;  there 
are  2*"  sets  in  the  domain  of  the  transformation  that  are 
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transformed  into  a  single  set  in  the  image.  For  all  cases 
the  absolute  value  of  the  Jacobian  of  transformation  is 
one.  The  joint  density  of  A  ,  ...,  then  becomes 

^  t-1 

f(yi» . ^t-i'®^  ^  2. 

i=l  k=l 

where  X2»  ...»  x^)  is  the  p.d.f.  of  a  multivariate 

normal  distribution,  and  for  every  i  one  combination  of  the 

A 

signs  is  valid.  The  distribution  of  A  then  naturally  re¬ 
quires  the  (t-1) -fold  integration  of  the  joint  density  of 
A  ,  Y^,...,  •  Although  every  parameter  in  the  ]oint 

density  is  known,  the  labor  required  to  arrive  at  an  ex¬ 
plicitly-defined  probability  distribution  makes  the  task 
very  difficult.  Therefore,  the  analysis  concerning  the 
dispersion  of  a  was  not  completed  because  of  the  intracti- 
bility. 


25 


IV.  MOMENTS  OF  & 

The  next  step  of  the  analysis,  then,  was  a  simulation 
of  the  inventory  system  in  order  to  gain  an  insight  concern¬ 
ing  the  distribution  of  ^  .  It  was  seen  that  the  reorder 
levels  were  fluctuating  around  the  theoretical  reorder  lev¬ 
els,  even  after  1000  observations.  This  was  related  to  the 

asymptotically  finite  variance  of  A  as  indicated  by  [Ref.?]. 

❖ 

Although  computer  simulations  implied  variance  of  A  does 


not  vanish,  the  moments  of  A  were  derived  in  order  that  an 
analysis  of  the  simulation  results  might  be  made.  The  mo- 

A 

ments  of  A  were  derived  for  large  values  of  t  and  are, 
therefore,  asymptotic. 


Cov(|e^_^I,  lej._j|)  =  E(l  ^  ^  "  ^<1  ®t-i^  *t-j^  ^ 

(4-2) 

Thus,  substituting  (3-3),  (3-4),  and  (4-2)  into  (4-1), 


Vih  =  2  “ 


2-Ot  I+I 


M. 


(1+1*) 


(t2  + 


t-2  t-1 


2“  2.  2, 

i=0  j  =  '!+l 


(4-3) 


This  last  expression  for  V( A)  requires  the  computation  of 
E (I  e^_^e^_j  1  )  .  The  following  derivations  concern 


E<l  n-iet-jl)- 


The  joint  distribution  of  indicated 

in  Section  III,  is  bivariate  normal.  The  parameters  of 
this  distribution  are 


Ul  =  E(ej._^)  =  0 


ri  -  o 


>■2  =  -  0 


=  a  .1  + 


'2  -  ■  l+|i 


and  asymptotic  p  is  given  by  (3-2)  .  For  notational  con¬ 
venience,  let 


X  =  e 


t-i 


Y  =  e 


t-j 


1 


i  ' 


f 

T 

i 


N 


t 

.4 


i 
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Then  the  joint  p.d.f.  is  given  [Ref.  4]  by 


f(x,y)  = 


_  a 
2 


2n  a  a 

X  y 


«  00  ^  ^  ^  00 

—  ®  <  y  <  00 


where 


q  = 


l-p  a 


X  y  (7, 


In  that  case 


00  00 


E(1xy1  )  =  j  1  I  xyl  f(x,y)dy  dx 


«oo  _ao 


U  CD 

=  [  i  x|  T-l  y  f(x,y)dy  +  y  f(x,y)dy')  dx 
_00  .00  ~  ^ 


O  O 


=  I  J  xy  f(x,y)dydx  -  J  J  xy  f(x,y)dydx  - 


,00  ,0  00  00 

J  J  xy  f(x,y)dydx  +  J  J  xy  f(x,y)dydx 


o  o 


However,  by  the  definition  of  correlation  coefficient. 


o  o 
0  0 


CD 


E(XY)  = 


xy  f(x,y)dydx  +  [  [xy  f(x,y)dydx  + 
00  -'-oo  •Jq 

j  J  xy  f(x,y)dydx  +  j  [  x  y  f(x,y)dydx 


_00  _00 


O  -® 


O  O 


=  a  a  p 
x  y 
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Thus, 


o  o 


E (I  xy| ) +  E(XY)  =2^1  ^  xy  f(x,y)dydx  +  J  xy  f(x,y) dydx  ) 

mOD  mOB  —  — 


Making  a  change  of  variable  in  the  first  integral,  let 


z  =  -x 


u  =  -y 


then 


o  o 


J  xyf(x,y)dydx  =  [  zu  f(-z,-u)dudz 


o  o 


Since  f(x,y)  =  f(-x,-y),  it  follows  that 


-O  .o 


j*  xyf(x,y)dydx  =  J  |  zu  f(z,  u)dudz 


.00  .ao 


o  o 


Thus 


00  00 
ft  p 

E(Ixy|)  +  E(XY)  =  4  J  xy  f(x,y)dydx 


(4-4) 


The  right-hand  side  of  (4-4),  without  the  coefficient  4,  is 

x^  (v-b) ^ _ 

•  "  ”  2(3^  “  2f7^^(l-p^) 

[  J  xyf(x,y)dydx  =  I  — -  j  - -  dy dx 

O  O  O  CTxv^2ff  CTy/l^ 

(4-5) 

where  b  is  the  conditional  mean  of  Y  given  X,  and  is  given 


by  [Ref.  4], 


b  -  ^  px 
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The  conditional  standard  deviation  of  Y  given  X  is  given  by 

t  A 

c  =  a  n/I-p 

y 

Let 


z  ,  then 


o 


_  -(Y-b)2 

iQ- 


dydx  = 


xe  ^  r  (cz+b) 

■'o  Jb 


e 


dz  dx 


Integrating  the  inner  term  first: 


(4-6) 


cz+b 
— e 

-b 

c 


c 


(1 


where  «I>(x)  is  the  cumulative  normal  distribution.  Sub¬ 
stituting  this  expression  into  (4-6), 


I 


2a 


xe 


o  a 

X 


I 


/  ^2.  2 

_  (y-b  ) 

“  2c^ 

ye _ 

o  c  y2ff 


2a 


dydx  = 


cx  e 


o 

2 


a  X  2»r 


2a 


2a 


bx  e 


o  a  ./2ir 

X 


dx  - 


'  D  9(-  —  )  e 
c 

o  a  v'2»r 

X 


dx 


2c 


dx  + 


(4-7) 


30 


T 


m 


Denote  the  terms  on  the  right-hand  side  as  1,  2,  and  3, 


respectively,  and  compute  each  one  separately  as  follows: 


1.  Substituting  values  of  c  and  b  into  1, 


2 

P  X 


cx  e 


2(i-p2)tyx^ 


y - J  20  ^a-P^) 

or.yi-p  X  e 


ff.Vi-P 


IF—  "xfi-o  ' 


a  a 

—  V  X 

2ir 


2.  Substituting  the  value  of  b  into  2, 


2 

.»  My  e  ^ 


“  -  ax  =  !4-  f'x^  e 

o  a  yiF  cr  ^v/2ir  •’o 

X  X 


a  a  p 
y  X 


3.  Substituting  the  values  of  c  and  b  into  3, 


I  Ht^) 


o  a^Vl-p 


v/l-p^'  ^  (Tj2n 


20  ‘ 

X  o 

X  e  V  - 

- — -  ^  X  p  = 

O  v^2ff  X 

X 


o  ^J2n  '  o 


2 

■  ^  2  -?  _  tl 

2  ,  e  2 

j_.  721'' 


where 


{4-8) 


0^1 -p'^ 


J 
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Let  3c/cr^  =  y  in  (4-8) ;  then  (4-8)  becomes 


a  a  p 

XV 

j2iP 


o 


dt  dy 


where 


Integrating  by  parts. 


Eq.  (4-8)  becomes 


K 


-y 


dt 


00  00 
n 


dt 


dy  . 


where 


The  first  term  in  the  above  equation  is  zero;  hence, 
(4-8)  becomes 


2w 


+  K 


r 


dt  dy 


(4-' 


Now  considering  the  second  term  in  (4-9) , 


(4-11 

This  result,  due  to  the  fact  that  y  changes  over  positive 
numbers  anct  p  la  between  -.5  and  zero,  makes  the  argvunent 
of  ♦  (•)  positive.  Hence  (4-10)  becomes 


®  ^dtdy  =  M 

o  ■  o 

Let 

t  =  V  y? 

y  = 

then 


where 


JT  dv  du 


2 


n=0 


Since  the  series  is  uniformly  convergent,  it  fallows  that 


00  00 
(*  ^-u 


2n+l 


M  = 


2  2  y2ff  ^  ^ ^ 


u 


n  +  i 


o  u 


n=0  (^/l-p  ) 


2-  2n+l  n!(2n+l) 


K-." 


l=£l 


2n+l 


n  -u  , 
u  e  du 


n=0 


/ - 9  2n+l  2n+l 

(^/l-p^) 


n: 


or. 


2n-H 


M 


n=0 


•  2n+l 


It  was  previously  shown  that  p  is  always  between  -.5  and  0 
Therefore, 


-  1  < 


_  _fi- 


<  1 


Hence  the  series  in  the  above  expression  is  the  expansion 


Thus  Eq.  (4-10)  becomes 
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/  ot 

But  (T^  and  are  both  equal  to  (T  ^  1  +  .  Hence, 


e(|xy|  )  = 


+  7^) 


-  P  Arctan  (y  ^  ^■))  (4-13) 


Therefore, 


M. 


t-2  t-1 


(1+P‘)  »r(l-j9) 


i=^  j=d.+l 


■J 

(4-14) 


where  ® ( I ®t-j®t-i ^ ^  given  by  (4-13). 

A 

It  is  clear  that  to  arrive  at  an  expression  for  V(A) 
in  the  form  of  (4-14)  many  tedious  and  difficult  calcula¬ 
tions  must  be  made.  In  fact,  in  order  to  arrive  at  an 

A 

explicit  value  for  V(A),  one  must  deal  with  an  arctangent 

2  . 

term.  However,  because  cr  is  common  to  each  term,  the 
ratio  V(A)/a  can  be  defined  in  terms  of  the  known  con¬ 
stants,  and  an  approximate  numerical  answer  can  be  found 
for  the  ratio  V(A)/(7  .  'rne  computer  program  in  Appendix 
B  computes  the  first  two  terms  in  (4-13)  exactly,  and  the 
third  term  approximately.  Briefly,  the  computation  of  the 
third  term  was  accomplished  as  follows; 

Eq.  (4-13)  is  a  bounded  number  for  a  given  a  and  j8., 

It  cannot  be  made  greater  than  some  finite  number  d,  what¬ 
ever  be  the  value  of  p  .  (Note  that  p  is  the  only  variable 
in  (4-13)  that  depends  on  i  and  j,  and  it  is  bounded.) 
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Therefore,  in  the  double  summation  of  (4-14),  for  some 
finite  value  of  i  and  j,  higher  order  terms  will  get  smaller 
than  a  specified  number  m  which,  in  this  case,  was  chosen 
as  10  The  program  was  written  so  that  the  higher  terms 

smaller  than  were  not  included  in  the  result,  and 

were  subsequently  dropped  from  the  summation.  The  upper 
bound  of  the  error  introduced  by  omitting  the  higher  terms 
from  the  summation  can  be  found  as  follows: 
e(Ixy|  )  is  bounded  from  above  by 


Since  the  summation  is  infinite, 

00  00 

Upper  Bound  of  Error  =  20?‘  lO”^®  ^ 

i=0  j=d+l 


First,  this  error  is  additive  and  hence  is  added  to 
the  result.  Second,  whatever  be  the  value  of  P  ,  it  is 
always  much  smaller  than  one.  Hence,  the  results  contained 
in  Table  1,  obtained  by  the  computer  program  for  various 
values  of  a  were  almost  exact. 
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SMOOTHING 

CONSTANT 

VARIANCE  OF 

STANDARD 

DEVIATION  OF  (a)/ff 

0.10 

0.0204 

0.1745, 

0.15 

0.0325 

0.2173 

0.20 

0.0461 

0.2553 

0.25 

0.0614 

0.2905 

0.30 

0.0785 

0.3239 

0.35 

0.0979 

0.3561 

0.40 

0.1196 

0.3876 

0.45 

0 . 1440 

0.4187 

0.50 

0.1716 

0.4495 

0.55 

0.2026 

0.4802 

0.80 

0.4267 

0.6341 

Table  1.  Variance  of  — and  Standard  Deviation  of  —  . 

<y 
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^  COMPARISON  OF  EXPONENTIAL  SMOOTHING 
AND  MLE  BY  COMPUTER  SIMULATION 

For  the  reason  explained  previously,  simulation  was 
utilized  to  explore  the  parts  of  the  model  where  theoreti¬ 
cal  work  was  very  difficult,  and  was  also  used  to  coopare 
the  exponential  smoothing  with  MLE.  The  reason  that  MLE 
was  chosen  as  the  alternative  procedure  was  because  it  has 
the  property  of  an  asymptotically  vanishing  variance.  As 
indicated  in  [Ref.  s],  the  maximum  likelihood  estimator, 

A 

£7  ,  of  C  is  given  by  the  formula 


A 

a 


I 


i=l 


(x.-x) 

n 


2 


where  x  is  the  sample  mean.  The  statistic  a  has  a  chi 
square  distribution  [Ref.  3];  in  fact. 


T{%) 


n^) 


and 


r(Y^) 


From  these  formulas. 


Lim  E{ct) 
n-“ 

Lim  ^„V(a) 

n-oo 


a 

0 
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The  computer  program  used  for  the  simulation  is  included  in 
Appendix  A.  The  simulation  procedure  and  ensuing  results 
are  as  follows: 

One  thousand  normally-distributed  random  numbers,  re¬ 
presenting  the  demand,  were  generated  by  a  pseudorandom 
number  generator.  Various  values  of  the  parameters,  4,  d, 
and  ill  were  used.  For  evex*y  triplet  of  the  parameters,  re¬ 
order  levels,  risk  levels  attained,  estimates  of  c.  and 
sample  mean,  and  sample  standard  deviation  of  d  and  of  the 
estimate  of  d  by  MAD  were  computed.  The  value  of  was  0.2 
for  all  cases.  Table  2  presents  the  results  achieved  by 
the  simulation.  As  Table  2  indicates,  the  theoretical 
values  of  the  standard  deviations  of  both  estimates  of  d 
were  in  complete  agreement  with  the  sample  standard  devia¬ 
tions  of  these  same  estimates.  The  variance  of  d  is  higher 

A 

than  that  of  d  in  every  case.  The  last  two  columns  give 
the  mean  square  error  (MSB)  of  d  and  d.  The  MSB  is  given 

A 

by  the  variance  plus  the  square  of  the  bias.  Although  d 
is  biased  for  1000  observations,  its  bias  can  be  neglected. 
Tnen  the  MSB  of  d  and  d  are  determined  by  their  variances. 
The  comparison  of  the  MSB  columns  are  quite  striking;  the 
MSB  for  d  reaches  values  as  high  as  100  times  that  for  d. 

By  noting  the  similarity  between  the  theoretical  values  and 
the  simulation  results,  the  effect  on  reorder  levels  was 
analyzed  as  follows: 

Computation  of  reorder  levels  was  made  by  utilizing 
Eq.  (3-3).  First,  the  theoretical  values  of  the  reorder 
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Table  2,  Seunple  Characteristics 


levels  were  computed  for  every  fixed  parameter  triplet  fi, 
a,  j|)  .  Then  the  reorder  levels,  found  by  smoothing  and 
MLE,  were  compared  with  the  theoretical  values.  This  com¬ 
parison  is  summarized  in  Table  3.  The  value  in  Table  4 
is  the  percentage  (rounded  off  to  the  nearest  integer)  of 
1000  reorder  levels  found  to  be  within  one  unit  of  the 
theoretical  reorder  level.  The  first  colvimn  under  each  of 
the  headings  and  lists  the  results  for  smooth¬ 

ed  estimates;  the  second  columi)  under  each  heading  lists 
MLE.  The  quantity  represents  the  number  of  units  about 
the  theoretical  reorder  level  within  which  50  per  cent  of 
the  computed  reorder  levels  were  found;  Pg^  is  a  similar 
computation  for  95  per  cent  of  the  reorder  levels.  Table 
3  clearly  illustrates  the  superiority  of  MLE  ccxnpared  to 
exponential  smoothing.  For  the  case  where  smoothing  was 
used  to  estimate  CT,  reorder  levels  fluctuated  around  the 
theoretical  reorder  level,  causing  either  an  excessive 
purchase,  or  imposing  more  risk  than  desired.  The  magni¬ 
tudes  of  these  fluctuations  are  clearly  evident  from  the 
table  . 

Table  4  indicates  the  risk  levels  attained  by  both 
methods  at  the  1000th  observation.  For  every  predetermined 
parameter  pair  n  and  cr,  five  different  risk  levels  were 
selected.  For  every  triplet,  the  theoretical  reorder 
levels  were  compared  with  the  computed  reorder  levels  of 
both  procedures.  The  actual  risk  levels  attained  were  then 
computed. 


42 


PARAMETER 

TRIPLETS 

Qi 

^0 

^5 

40,4, .01 

27 

95 

2 

0 

6 

0 

40, 4,. 11 

41 

94 

2 

0 

4 

0 

40,4, .50 

49 

92 

1 

0 

3 

1 

50, 5, .01 

22 

96 

3 

0 

8 

0 

50, 5,. 11 

32 

98 

2 

0 

5 

0 

50, 5, .50 

46 

99 

1 

0 

3 

0 

100.10, .01 

12 

90 

5 

0 

14 

1 

100, 10,. 11 

18 

94 

4 

0 

11 

0 

100,10, .50 

26 

94 

2 

0 

6 

0 

200,20, .01 

5 

71 

10 

1 

29 

6 

200,20, .11 

9 

82 

6 

1 

18 

3 

200,20, .50 

10 

84 

5 

1 

14 

2 

300,30, .01 

4 

81 

12 

1 

37 

6 

300,30, .11 

7 

74 

8 

1 

24 

3 

300,30, .50 

8 

51 

7 

1 

18 

3 

400,40, .01 

2 

15 

20 

3 

53 

10 

400,40, .11 

4 

32 

13 

2 

37 

6 

400,40, .50 

5 

75 

10 

1 

27 

3 

500,50, .01 

3 

1 

24 

9 

58 

11 

500,50, .11 

3 

4 

15 

4 

38 

6 

500, 50,. 50 

4 

73 

11 

1 

29 

4 

600,60, .01 

1 

12 

27 

3 

77 

17 

600,60, .11 

3 

10 

18 

3 

51 

5 

600,60, .50 

3 

1 

14 

4 

42 

12 

700,70, .01 

1 

33 

32 

3 

93 

16 

700,70, .11 

3 

35 

20 

2 

63 

10 

700,70, .50 

4 

16 

15 

2 

45 

15 

800, 80,. 01 

1 

0 

37 

12 

98 

27 

800,80, .11 

2 

0 

25 

8 

64 

18 

800,80, .50 

3 

32 

18 

2 

54 

9 

Table  3.  Percentiles  for  Reorder 
Levels. 


43 


PARAMETER 

TRIPLETS 


PARAMETER 

TRIPLETS 

^50 

^95 

900,90, .01 

1 

30 

42 

2 

125 

B 

900,90, .11 

2 

37 

29 

2 

83 

B 

1000,100, .01 

1 

25 

45 

3 

126 

12 

1000, 100, .11 

1 

8 

29 

4 

87 

1000, 100, .50 

1 

6 

25 

2 

68 

Table  3 . (Continued)  -  Percentiles  for  Reorder 

Levels. 


1 

! 

I 


As  Table  4  indicates,  the  risk  levels  attained  by  MLE 
are  found  to  very  close  to  the  desired  rlak  levels. 
However,  the  rlak  levels  of  the  sinioothlng  method  were 
found  to  fluctuate  over  a  wide  range. 

These  results,  gained  by  the  simulation,  indicated 
that  the  MLE  technique  la  to  be  preferred  over  the  tech¬ 
nique  of  exponential  smoothing.  However,  due  to  the  lack 
of  theoretical  analysis  concerning  the  distribution  of 
these  results  are  not  considered  conclusive.  However,  this 
analysis  does  indicate  that  whenever  a  constant  mean  de¬ 
mand  can  be  justified,  the  MLE  technique  is  to  be  prefer¬ 
red  to  the  exponential  smoothing  technique. 
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APPENDIX  A 


SIf'ULATION 


IM  the  Ff^LLOWING  PROGRAM  TWO  MFthqDS  OF  FSTIMSTION  OF 
MEAN  ANO  STANDARD  DEVIATION  0^  A  STATIONAPY  ORGCESS 
APF  CONTOASTFO  BY  USING  SIMULATION  TECHNIQUES 
METHODS  ARE: 

1. THF  exponentially  SMOOTHED  DATA  AS  THE  ESTIMATF  OF 
MEAN  AND  the  EXPONENTIALLY  SMOPTHEO  ABSOIUTE  VALUES 
OF  EPROPS  AS  THE  ESTIMATE  OF  STANDARD  DEVIATION  TIMES 
A  KNOWN  CONSTANT 

2.  MAXIMUM  LIKELIHOOD  ESTIMATES 

ALSO  IN  THE  PROGRAM  THE  EFFECTS  OF  BOTH  ESTIMATION 
PROCEDUPFS  ON  THE  RISK  AND  OFOROER  LEVELS  ARE  SIMULA¬ 
TED  BY  S/MPLE  SI7E  1000 


DIMENSION  0K(6)  ,T(6I  ,XL(6)  ,XP(6 ) iXM(6)XN(6) ,XK(6  » 
IXK  P0Q,5)  ,X2<  1000.51  ,G1  (1005)  .G2  (10  05  )  ,  ?  ( 10  ) ,  TX  ( 6 ) 
2TY(6)  ,0ISKX(6)  ,PISKY(6) 


RISK  LEVELS 


PK( 1)=0.01 
eK( 2)=C.C5 
PK( 3)*C.ll 
eK( 4)=c.25 
0K( 5)=C.fC 

RFAn(5,131l  (T(I),I=1,E) 
131  E0PMAT(5F4.2) 

S=10.Q 

XJ=C.C 

Y*URN(C) 

no  021  MT=1,16 
XJ=XJ  +  1.C 
no  1  j=i,5 
DO  1  1=1,800 
Xl(  I,  JI--C.0 
1  X2(I,J)=C.O 
?L=S 


smoothing  constant  alpha  is  3,2 


BFTA=0.P 

WRITF(ft,87S)  S.STPFX 

PTC  COPMA  T(  *1  'tA5,  'MEAM^*  ,P6.1  ,nx,  '^.pcvsi  ,F6.1) 
OELT^  =C.7S8*SORT(2.C/(1.0+BETA)  )*FTDFX 


rOMOliTATION  OF  THE  THEORETICAL  PEOPDER  LEVELS 


no  11  1=1,5 

11  7( I )=  7L*ST0EX*T(I  ) 

WRTTE(6,15C)  (7(11,1=1,5) 

15C  POPMAT( //TlO.'THFnRE'^ICAL  REORDFP  LEVELS ',  SF  10.  3/7  ) 
v=i,r. 

SUM=7L 

TnTxo,o 

R  1=0,0 

R2=0,0 

B3=C.C 

PA=C.C 


GFMEPATIGN  OF  NORMALLY  DISTIBUTED  NUMBERS  WHICH  HAVF 
MEAN  S  AND  'JTANDARD  DEVIATION  STDEX 
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DO  21  1=1, loco 

Cfllt  NORMAL  IC.O.STDFXfX) 

Y=ZL+X 


c  method  1 


DELTA=0.2*ABS(S-Y)+C.8*DELTA 

S  =  C.2*Y+C.'i  +  S 

RDEV*l.iePTDELT6 

V=VT1.0 

S»1M=<L'M+Y 


C  METHOD  2 


X8AR=SUM/V 

TOT=TOT+! XBAR-Y)**2 

STnEV.’:SOPT(  TCT/(  V-l.C)  I 

r,l(  I)  =  SPEV 

G2(  I  )  =  ?TOEV 

on  20  J=l,5 

XL( J»=S+S0GVTT( J» 

XP(.I)=XRAR  +  STDEV*T{  J) 

XK( J)=Xl( J)-Z( JJ 

TX( J)=(XL( J)-ZL)/STDEX 


C  RISKS  ATTAINED  BY  BOTH  METHODS 


RTSKX( J)=C. 5-0. 5TERF(TX(J) *0.707) 

TY(  J)  =  (XP(  Jl-ZDZSTDFX 

PISKV( J)=C. 5-0. 5*EPF(TY( J| *0.707) 

XM( J)sXF( J)-Z( J) 

N1  =  XK(  J) 

N2*XM(  J) 

I!=(Nl.LT.O)  Nl=-Ml 
IF(N2.LT.O)  N2=-N2 
XI{NU1,J)*XUN1  +  1,J)*1.0 

20  X2» N2T1.J)=X2(N2+1 ,J)+1.0 
NS=I 

IF(  NS.LF.ePa)  GO  TC;  21 

WR !TE(6,1CI)  ( XL(M) .M=l  ,5)  ,(XP(M) ,M=1,5 ) 

ICl  FORMAT! /T13,5F7,1 ,10X,5F7.l) 

WRITE (6, 147)  (RISKX(M) ,M=l -5) .(RISKYlM) ,M*l,5) 
1A7  FORMAT! /T2, 'ACTUAL  R  T  SK  »  ,5  F7,3  ,  lOX  ,  5  F7 , 3  ) 

21  CONTINUF 


f  COMPUTATION  OF  THE  OBSERVED  RRFQUENCIER  WITHIN  THE  I 
C  NEIGHBORHOOD  OF  ACTUAL  REORDER  LEVEL  WHERE  1=1, 2, 3... 


DO  26  J=l,5 
WRITF!6,105)  PK!J) 

1C*  FORMAT!//T45,'PISK=',F5.2) 

A  1=0.0 
A  2=C.C 
1  =  1 

25  IF! Al.GF.  1.0)  GO  TO  24 
A l=Xl! I ,J) /ICOO.O+Al 
24  IF!A2.GE.1.0)  GO  TO  35 
A2=X2! T,J)/1CC0.0+A2 
35  WRITE!6,112)  I,A1,A2 

112  ROPMAT! TiO, 'PERCENTILE  BETWEEN  thEO.  REORDER  LEVEL 
1  ♦AND-',I4,2F10.4) 

1  =  1  +  1 
XA=A  1+A2 
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IF(X4.LF.l.gS6)  GOTO  25 
26  CONTINUF 


C  SAMPLE  *'EAN  AND  STANDARD  DEVIATION  OF  BOTH  ESTIMATES 

C  OF  POPULATION  STANDARD  DEVIATION 


DO  I"'  r*i.icoo 

ei=r,i(  I  )+Bi 

27  B2=02+G2«II 
B1=RI/1C00.0 
B2=R2/ICC0.0 
DO  2R  1=1,1000 
B3=B3+(G1(  I I-B1)**2 

28  8A=84*<B2-G2m  l♦♦2 
B3=SnRTfR3/lOOO.O) 

B4=SORT(B4/1COO.O) 

WRITF(6,200) 

200  PORMAK //T45, 'SAMPLE  MEAN  AND  S.  DEV  OF  ESTIMATED  S. 
IDEV  I 

WRITE(6,20l)  B1,83,B2,R4 

2C1  pOPM/T(//T4*>,'(MADI  '.SX.'MEAN*'  ,F7. 2  , »  S  DEV=  • »  F7 , 2»  1 5X  , 
l'( S2) • ,fX, •MFAN=»,F7.2,  'SDEVx* ,F7,2I 
S=10C.C 

IF(XJ-5.C)  72l,721tT22 

721  S=XJ*1C.0 
GO  TO  821 

722  S  =  S*(XJ-5.0) 

821  CONTINUE 

STOP 

END 


SUBPnUTTNe  normal  (EX,ST0X,X» 

SUM=O.C 

DO  5  1*1,26 

R=URN( 1) 

SUM=SUM+P 

SUM=SUM/26.0 

X=STDX*17.614(SUM-0.5|+EX 

return 

END 


APPENDIX  B 


numfotc/m  computation  of  vaomncf  of  a 


IMPLICIT  (A-H,C-Z) 

RFT1  =  .0 

nn  20  K=i,ii 


SM^^THT^o  cgnft^mt 


ALPHA  =  1.-BETA 
A'5UM  =  0.C 

wqiTC(6tllO)  At  PH A 

li:  FGRMAT(//T1C,  “^MOCTHIKO,  CONSTANT  =  •  ,  . 3  ) 


apcpc  XIM.A  T TPN  rp  thF  INFINTTP  '1)M 


no  1  M=l,3G3 
I=M-1 

J  =  M 

T=/i  .0/1  3.  1,0+flFTAn 

P  =  Al  PHA/(  I.+BE-^A  ) 
n=BCTA**(  UJ) 

C=BCT4^*(  _,|.i  ) 

4  =C  /P  E  re 
3  B=A»eE‘^A 

RHn  =  I  l.+BETA  )/-- .  ♦(  1.  +  F*(1.  u-A**2  +  (  1  .  ^BI  ♦«2)  +  (  l.~ ALPHA* 
1AI**2-1. 

7N=nSQR''(  l.-RHn**2t 
7C  =  _rhG/ZN 

<!UM=T*(  7N-PH0*nATAN(  ZP))  *0 

ASU^*=ASL'M+'I|W 

IP(  SUM-1, OP-10)  Itl,? 

2  A=A*BETA 
n=D*BETA 
GO  Tn  3 

1  continue 

WPITr(6,iC2)  SUM 

1C?  fOPMAT(T1C,‘LAST  SUM=SF20.12) 

X=ALCHA/(  l.  +  RFTA)  ♦*2*2.282 /3.1'-1 
Z=P.C /3. 141*RFTA /( 1.0+RFTA) **2 


VAR lA  NCF  0^  A 


VAP  =X-.7+2.  G*(  1, 0-RFTA  )  ♦*2*AFiJM 


ST/MHAPn  P'=VIA''ICN  PP  a 


cnp  v=G  SOP  3.1^1  ♦(  l.  +  R'^TAI  ♦VAR)  /2.C 

wR ! te ( 4 , i c I )  VAP, snt V 

ICl  FORMA  TITIO. 'VARIANCE  PF  MAO  •  ,F1 .7 , 5X  ,  '  S'"  ANPART  CEVIA 
ITTOM  •tFlE.T) 
flE^A  =  PFTA-.05 
IF(K.FO.IC)  BETA=.2 
20  CONTI NUF 
cTnp 
FNP 
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